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The boundary integral method is extended to derive closed 
integro-differential equations applicable to computation of the 
shape and propagation speed of a steadily moving spot and 
to the analysis of dynamic instabilities in the sharp bound- 
ary limit. Expansion of the boundary integral near the locus 
of traveling instability in a standard reaction-diffusion model 
proves that the bifurcation is supercritical whenever the spot 
is stable to splitting, so that propagating spots can be stabi- 
lized without introducing additional long-range variables. 



Localized structures in non-equlibrium systems (dis- 
sipative solitons) have been studied both in experi- 
ments and computations in various applications, includ- 
ing chemical patterns in solutions Q and on surfaces , 
gas discharges [[| and nonlinear optics . The interest to 
dynamic solitary structures, in particular, in optical Q] 
and gas discharge systems ^ has been recently driven by 
their possible role in information transmission and pro- 
cessing. 

A variety of observed phenomena can be reproduced 
qualitatively with the help of simple reaction-diffusion 
models with separated scales ]6|-|l0| . Extended models of 
this type included nonlocal interactions due to gas trans- 
port §,0, Marangoni flow Q or optical feedback ||,|f3|. 
A great advantage of scale separation is a possibility 
to construct analytically strongly nonlinear structures in 
the sharp interface limit. An alternative approach based 
on Ginzburg-Landau models supplemented by quintic 
and/or fourth-order differential (Swift-Hohcnberg) terms 
|l4| have to rely on numerics in more than one dimension. 

Dynamical solitary structures are most interesting 
from the point of view of both theory and potential ap- 
plications. Existence of traveling spots in sharp-interface 
models is indicated by translational instability of a sta- 
tionary spot [pd| . This instability is a manifestation of 
a general phenomenon of parity breaking (fsing-Bloch) 
bifurcation [fT5|Jl6l which takes a single stable front into 
a pair of counter-propagating fronts forming the front 
and the back of a traveling pulse. Numerical simula- 
tions, however, failed to produce stable traveling spots 
in the basic activator-inhibitor model, and the tendency 
of moving spots to spread out laterally had to be sup- 
pressed either by global interaction in a finite region |Tl| ] 
or by adding an extra inhibitor with specially designed 
properties Ji7[. 



The dynamical problem is difficult for theoretical 
study, since a moving spot loses its circular shape, and a 
free-boundary problem is formidable even for simplest ki- 
netic models. Numerical simulation is also problematic, 
due to the need to use fine grid to catch sharp gradients 
of the activator; therefore actual computations were car- 
ried out for moderate scale ratios. A large amount of 
numerical data, such as the inhibitor field far from the 
spot contour, is superfluous. This could be overcome if it 
was possible to reduce the PDE solution to local dynam- 
ics of a sharp boundary. Unfortunately, a purely local 
equation of front motion |jl6f is applicable only when the 
curvature far exceeds the diffusion scale of the long-range 
variable, whereas a spot typically suffers splitting insta- 
bility |^| before growing so large. On the other hand, 
the nonlocal boundary integral method (f8) is applica- 
ble only when the inhibitor dynamics is fast compared 
to the characteristic propagation scale of a front motion, 
i.e. under conditions when no dynamic instabilities arise 
and traveling spots do not exist. 

It is the aim of this Letter, to extend the nonlocal 
boundary integral method to dynamical problems, and 
to find out with its help conditions of supercritical bifur- 
cation for steadily moving spots. We consider the stan- 
dard FitzHugh-Nagumo model including two variables - 
a short-range activator u and a long-range inhibitor v. 



e 2 Tu t = e 2 V 2 u + V'{u) - ev, 
Vt = V 2 t> — v — v + /iu. 



(1) 
(2) 



Here V(u) is a symmetric double- well potential with min- 
ima at u = ±1; e < 1 is a scale ratio, and other parame- 
ters are scaled in such a way that the effects of bias and 
curvature on the motion of the front separating the up- 
and down states of the short-range variable are of the 
same order of magnitude. The local normal velocity of 
the front is 



^(bv- /c) + 0(e), 



(3) 



where n is curvature and b is a numerical factor depen- 
dent on the form of V(u); for example, b = 3/V2 for the 
quartic potential V(u) — — 4(1 — u 2 ) 2 . By definition, the 
velocity is positive when the down-state u < advances. 

In the sharp boundary approximation valid at e <C I , 
a closed equation of motion for a solitary spot propagat- 
ing with a constant speed can be written by expressing 
the local curvature in Eq. (^) with the help of a suit- 
able parametrization of the spot boundary, and resolving 
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Eq. (J2J) rewritten in a coordinate frame propagating with 
a speed c (as yet unknown). It is convenient to shift the 
long-range variable v — w—f+fi, so that w(po) = when 
the up-state u = 1 — 0(e) prevails at infinity. The sta- 
tionary equation of w in the coordinate frame translating 
with the speed c is 



c • Vw + V 2 w - w = 2/j.H, 



(4) 



where, neglecting 0(e) corrections, H = 1 inside and 
H = outside the spot. The solution can be presented 
in the form of an integral over the spot area S: 



t»(x) = -£ f g(x-t)d 2 z, 



(5) 



where the kernel Q contains a modified Bessel function 
K : 



Q{v) = ^-e-^K [\v\^l + ^ 



(6) 



This integral can be transformed into a contour integral 
with the help of the Gauss theorem. To avoid divergent 
expressions, the contour should exclude the point x — 
Clearly, excluding an infinitesimal circle around this 
point does not affect the integral (||), since the kernel 
(||) is only logarithmically divergent. Replacing Q(r) = 
V 2 Cf (r) + c • V<?(r) (r ^ 0), we transform the integral in 
Eq. (g) as 

- J 0(x - d)d 2 d = J • H(x - Z)d 2 £ 

= j n(s)-K(x-$(s))ds, (7) 

where H(r) = V<?(r) + cQ(r) and n is the normal to the 
contour V . The vector Green's function H corresponding 
to the kernel in Eq. is computed as 



H(r) 



(8) 



When x is a boundary point, V consists of the spot 
boundary T cut at this point and closed by an infinites- 
imally small semicircle about x. The integral over the 
semicircle equals to n. Defining the external normal to 
r as the tangent t = x'(s) rotated clockwise by 7r/2, 
the required value of the long-range variable on the spot 
boundary (parametrized by the arc length s or a) is ex- 
pressed, using the 2D cross product x , as 

v(s) = I H(x(s) - x(cr)) x x'(a)dcr. (9) 

71" Jr 

To obtain a closed integral equation of a steadily mov- 
ing spot, it remains to define a shift of parametrization 



accompanying shape-preserving translation. Recall that 
Eq. (||) determines the propagation velocity c„ along 
the normal to the boundary. In addition, one can in- 
troduce arbitrary tangential velocity Ct which has no 
physical meaning but might be necessary to account for 
the fact that each "material point" on a translated con- 
tour is, generally, mapped onto a point with a differ- 
ent parametrization even when the shape remains un- 
changed. The tangential velocity can be defined by re- 
quiring that each material point be translated strictly 
parallel to the direction of motion, i.e. c„n + Cft = c. 
Taking the cross product with c yields c t = c„(c x t) 
/(c • t). Then eliminating q gives the normal velocity 
Cn = c x t necessary for translating the contour along 
the x axis with the velocity c. Using this in Eq. (|^) 
yields the condition of stationary propagation 



c x x'(s) = t 1 [bv(s) — k(s)]. 



(10) 



The form and the propagation speed of a slowly moving 
and weakly distorted circular contour can be obtained by 
expanding Eq. ( |l0| ) in c = |c| near the point of traveling 
bifurcation t = tq, which is also determined in the course 
of the expansion. For a circular contour with a radius a, 
Eq. (0) takes the form 



v(4>) = — v + 



fi.a 



^— ^ca (cos (p— cos ip) 



\ccosyK ( (2aJ 1 + 7-c 2 sin i|</> - tp 



■sin±|</>- <p\Jl + 



Ic 2 x 



K x 2a\ 1 



3 c 2 sin || 



<P\ 



dp, 



(11) 



where <f> or ip is the polar angle counted from the di- 
rection of motion. The angular integrals that appear in 
the successive terms of the expansion are evaluated itera- 
tively, starting from $0(0) = 7T ^o( a )^o( a ) an d using the 
relations 



sin 2fe+1 ^X 1 (2asin^)d^-i^, 
2 V 2/ 2 da ' 

1 dia^k^) 



$ fc (a) = j Bin 2 * I ^0(20 sin I 



2a da 



Effect of small boundary distortions on v can be com- 
puted directly with the help of Eq. (||) , where the integra- 
tion should be carried out only over a small area swept 
by the displaced spot boundary. This approach is most 
useful for stability analysis with respect to small pertur- 
bations of a known static shape, and is easier than using 
the expansion of Eq. (^J) with a perturbed boundary. For 
a circular spot, we expand the perturbations of both v 
and p in the Fourier series 

p(<f>, t) — p(4>) — a = c"a„e A "* cos n<fi, 

n>2 
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v(<j>,t) — ^ v n e Xnt cosncj). 

n>2 

The curvature is expressed as 
P 2 - 2 pI ~ PP<t><t> 



(12) 



<<t>) = 



(p 2 + pD 3/2 

a- 1 + 3(c/a) 2 a 2 e X2t cos 2(f) + 0(c 3 ). (13) 



Since the displaced point should remain on the bound- 
ary, the distortion p(tp) should be compensated by rigid 
displacement of the spot by an increment p{4>) when v(<fi) 
is computed (see the inset in Fig. [l]). The resulting equa- 
tion for eigenvalues A„ following from Eq. (|^) is 



n 2 — 1 Aabp 



cos ncj) d<f> x 
\p{tp) - p{(j>) cos(</> - tf>)]e- ? ca[cos ~ cos v) x 
K Q [2aJl + \„ + ±c 2 sin §|0-y>|) dtp. (14) 
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FIG. 1. The bifurcation diagram for stationary spots at 
T = 1. C - existence boundary, S - locus of splitting instabil- 
ity. The stability region is bounded by the locus of breathing 
instability B, branching off at the point of double zero eigen- 
value D, and the locus of traveling instability T. Inset: a 
circular spot distorted by second and third harmonics with 
amplitudes proportional to c™ . The shape is characteristic to 
a spot propagating to the right, and the amplitudes are chosen 
in such a way that the curvature on the back side vanishes. 
The center of the gray circle is shifted from the black to the 
gray spot to compensate the distortion at <f> = 0, so that the 
integral is taken over the area between the black contour and 
the gray circle when the effect of small distortions on the v 
field at this point is computed. 

Using the constant zero-order term in the expansion 
of Eq. ( |TT|) together with k = a -1 in Eq.(|3|) yields the 
stationarity condition 

v = -(ba)- 1 +(j,a[K 1 (a)I (a)-K (a)h(a)}. (15) 

A stationary solution stable against collapse or uniform 
swelling exists in the region in the parametric plane p, v 



(Fig. 0) bounded by the cusped curve C and the axis 
v — 0, p > 2/6. This curve is drawn as a parametric plot 
with v{a) given by Eq. ( |l5| ) and p(a) by Eq. (|t|) with 
c, n and Ao set to zero (or, equivalently, by the condition 
Fq(o) = 0, where F (a) is the right-hand side of Eq. Ji~5|). 

The first-order term in the expansion of Eq. ( pT| ) is pro- 
portional to cos <j), and should compensate at the travel- 
ing bifurcation point the left-hand side of Eq. (|l0]). This 
yields the bifurcation condition 

r - bpa[a(h(a)K Q {a) - I (a)Jfi(a)) + 2I 1 (a)Ki(a)}, 

(16) 

which coincides with the known result obtained by other 
means O]. The curve T in Fig. [l] shows the traveling 
instability threshold for To — 1. The static spot is un- 
stable below this curve; the locus shifts up (to smaller 
radii) as r decreases, and exits the existence domain at 
t < 1/4. At r > 1, the dominant instability at large radii 
is a static splitting instability. Its locus, determined by 
Eq. (O) with n = 2 and c = A2 = 0, is the curve S in 

Fig ET 

Another possible dynamic instability is breathing in- 
stability mHll. Its locus is s iven b y E T & witn 

c = n — and Ao = iui. The frequency to as a function 
of the spot radius a is computed by solving the equation 
tlo = a~ 2 Im F(a,u})/~Re F(a,ui), where F(a,u) is the 
right-hand side of Eq. (UJ) computed as 



F(a,Lu) = 2pa l\ {a\f\ + iLuj K\ [a\/\ + iiSj 

- J (aVl + iu) K (a^Jl + iiSj . (17) 

The curve B in Fig. |l| shows the bifurcation locus at 
t = 1. The instability region retreats to small radii 
(large v) at large t and spreads downwards as r de- 
creases. The balloon of stable solutions disappears al- 
together at r < 0.5 after the tips of both dynamic loci 
meet on the existence boundary. 

In the second order, Eq. ([ll]) yields a constant term 

„< a .°) = ~^a 2 [a(h (a)K (a) - 7 (a)#i (a) ) + h {a)Kx (a)] 

(18) 

and a dipole term -y( 2 ' 2 ) = q( 2 - 2 ) cos 2(f), where 

g(2,2) = l^W^K^a) - h(a)K Q (a)) 

-3/i(a)Jfi(a) + 2I 2 {a)K 2 {a)]. (19) 

The constant term is positive and causes contraction of 
the average radius of the moving spot by an increment 
a = -aW-V . 

The second-order dipolar term in the right-hand side of 
Eq. (|l0|), u( 2 > 2 ) = (f 2 ^a 2 cos 2(f>, as well as the third-order 
first harmonic term, u^ 3 - 1 ) = ^ 3,1 )a 2 cos(f>, needed for the 
solvability condition to follow, are read from Eq. (nj) 
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with n = 2 and A2 
order in c: 



0, respectively, in zero and first 



^ = 
5(3,1) = 



-3a~ 2 + 26/*[Ji(o)Jiri(a) - h{a)K 2 {a)l (20) 
bfj,a 2 Ix(a)Ki(a). (21) 



The coefficient g^ 2 ' 2 ^ vanishes at the splitting instability 
threshold (curve S in Fig. |l]) , and must be negative when 
the circular spot is stable. Consequently, the distortion 
amplitude is a 2 = — q( 2 > 2 )/q( 2 > 2 ) < 0, so that the dipole 
term causes contraction of the moving spot in the direc- 
tion of motion and expansion in the normal direction. 

Continuing the expansion to the third order, we com- 
pute the first harmonic term contributing to the solvabil- 
ity condition. The latter has the form fc — fee 3 , where 
t = t — To and the coefficient k determining the character 
of the bifurcation is computed as 



k = 



bn (V 3 ' 1 ) - 7£(o)oV a '°> - g^V^W 2 ^) ■ (22) 



The first term is the coefficient at the first harmonic in 
the third order of the expansion of Eq. ([TT|). The sec- 
ond term takes into account the second-order radius cor- 
rection to the first-order first harmonic term. The last 
term gives the effect of dipolar shape distortion; it be- 
comes dominant when the locus of splitting instability 
is approached. Stable traveling solution should be ob- 
served beyond the traveling instability threshold, i.e. at 
t < 0; hence, the condition of supercritical bifurcation 
is k < 0. The numerical check of the symbolically com- 
puted expression shows that the traveling bifurcation is 
always supercritical when the spot is stable to splitting. 
The traveling solution bifurcating supercritically must be 
stable, at least close to the bifurcation point where it in- 
herits stability of the stationary spot to other kinds of 
perturbations. 

The third harmonic term that appears in the third or- 
der of the expansion delineates, together with the second- 
order dipolar term, the characteristic shape of a translat- 
ing spot, pointed in the direction of motion and spread 
sidewise, as in the inset in Fig. 51 w hich has been also ob- 
served in numerical simulations |17| . Beyond the range of 
the bifurcation expansion, the shape, as well as the prop- 
agation speed can be determined by solving numerically 
Eq. (^) with v(s) given by Eq. (@) and curvature com- 
puted using the fully nonlinear expression in Eq. (|l3|). Al- 
though the boundary integral method reduces a PDE to 
a ID intcgro-differential equation, the equation is rather 
difficult. Iterative numerical solution (20) tends to break 
down rather close to the bifurcation point, as soon as the 
shape distortion becomes strong enough to flatten the 
spot at the back side. Since the boundary integral equa- 
tion is non-evolutionary, there is no way to distinguish 
between a purely numerical failure of convergence and a 
physical instability that would lead to lateral spreading 
observed in PDE simulations Jn]] . 



The above bifurcation expansion proves that a stable 
traveling solution does exist in the basic model (pi), (Q) 
in the sharp boundary limit. The result is applicable 
at 1 ^> c ^> y/e. It can be extended straightforwardly to 
models with more than one long-range variable, provided 
all long-range equations are linear. Stable traveling spot 
solutions should be, indeed, more robust in an extended 
model where they have been obtained in PDE simula- 
tions ]l7[ |, whereas in the basic model they require fine 
parametric tuning aided by the analytical theory. 
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